# #######################################################
# Simulation based on Twitter Mobilization Experiments 
# #######################################################
rm(list=ls())
library("estimatr") # ver 0.22.0
library("Formula")  # ver 1.2-3
load(file = "../data/Adj_sim.rdata")
source("estimator.R")
source("Twitter/simulation_setup.R")

# Generating Unobserved Network 
uProb <-  function(overlap, prob_g, prob_u_mean = 0.0075, n_unit){
  prob_u1 <- overlap*(prob_u_mean/prob_g)
  prob_u0 <- (1-overlap)*(prob_u_mean/(1-prob_g))
  uW <- matrix(NA, nrow = n_unit, ncol = n_unit)
  for(i in 1:nrow(uW)){
    prob_use <-  rep(prob_u1, times = n_unit)
    prob_use[W_use[i,] ==  0] <- prob_u0
    uW[i, ] <- rbinom(n = n_unit, size = 1, prob = prob_use)
  }
  return(uW)
}

# ##########
# DGP
# ##########
{
  
  sim_size <-  100
  k_sim_size <- 20
  est_G_l_f  <- est_T_l_f  <-  coef_mat_l_f <- list()
  for(h in 1:2){ # different outcome DGP
    cat(paste("h == ", h, ",", sep = ""))
    est_G_l  <- est_T_l  <-  coef_mat_l <- list()
    model  <- model_sim[h]
    for(z in 1:4){ # different overlap 
      cat(paste("z == ", z, ", ", sep = ""))
      est_T_all <- c()
      est_G_all <- c()
      coef_mat_all <- matrix(NA, nrow = 0, ncol =  4)
      
      for(k in 1:k_sim_size){
        cat(paste("k == ", k, ",", sep = ""))
        set.seed(10000*h + 1000*z + k*100)
        uW <- uProb(overlap = overlap_use[z], 
                    prob_g = prob_g,  prob_u_mean = 0.0075, n_unit = n_unit)
        uW_sum <- apply(uW, 1, sum)
        est_T <- c()
        est_G <- c()
        coef_mat <- matrix(NA, nrow = sim_size, ncol =  4)
        for(sim in 1:sim_size){
          set.seed(10000*h + 1000*z + k*100 + sim)
          treatment <- rbinom(n = n_unit, size = 1, prob = 0.5)
          G <- c((W_use %*% c(treatment))/W_sum)
          U <- c((uW %*% c(treatment))/uW_sum)
          # outcome 
          if(model == "linear"){
            Y <- intercept_l + coef_T_l*treatment + true_G*G + coef_U_l*U + 
              rnorm(n = n_unit, sd = sd_norm)
          }
          if(model == "interactive"){
            Y <- intercept + coef_T1*treatment + true_G*G + coef_T2*treatment*U + 
              coef_U1*U + coef_U2*(G*U) + rnorm(n = n_unit, sd = sd_norm)
          }
          
          # estimation 
          weight_use <- weight_IPW(treatment, G, W_sum, prob =  0.5)
          weight_use[weight_use > 10] <-  10
          lm_b <- lm_robust(Y ~ treatment*G, weights = weight_use)
          est_G[sim] <-  coef(lm_b)["G"]
          est_T[sim] <-  coef(lm_b)["treatment"]
          coef_mat[sim, 1:4] <- coef(lm_b)
        }
        est_G_all <- c(est_G_all, est_G)
        est_T_all <- c(est_T_all, est_T)
        coef_mat_all <- rbind(coef_mat_all, coef_mat)
      }
      est_G_l[[z]] <- est_G_all
      est_T_l[[z]] <- est_T_all
      coef_mat_l[[z]] <- coef_mat_all
    }
    est_G_l_f[[h]] <- est_G_l
    est_T_l_f[[h]] <- est_T_l
    coef_mat_l_f[[h]] <- coef_mat_l
  }
  output <- list("est_G_l_f" = est_G_l_f, 
                 "est_T_l_f" = est_T_l_f, 
                 "coef_mat_l_f" = coef_mat_l_f)
  save(output, file  = "../results/twitter_sim_output.rdata")
}
